# Import Python libraries
import pandas as pd
import numpy as np
import datetime as dt
import matplotlib.pyplot as plt
import requests
from bs4 import BeautifulSoup
import geopandas as gpd
from shapely.geometry import Point
import warnings
warnings.filterwarnings("ignore")

#Load raw data
emissions = pd.read_csv("co2_emissions.csv")
gdp = pd.read_csv("world_bank_gdppc_ppp.csv", skiprows=4) #skiprows to remove summary rows
hdi = pd.read_csv("undp_hdi.csv", usecols=range(37)) #usecols to remove hdi component columns
energy = pd.read_csv("energy-intensity.csv")
policies = pd.read_csv("climate_policies.csv")

# Rename and drop irrelavent columns
emissions['ISO'] = emissions['ISO 3166-1 alpha-3']
gdp['ISO'] = gdp['Country Code']
hdi['ISO'] = hdi['iso3']
emissions = emissions.drop(columns=['ISO 3166-1 alpha-3'])
gdp = gdp.drop(columns=['Country Code', 'Indicator Name', 'Indicator Code'])
hdi = hdi.drop(columns=['iso3', 'region', 'hdicode', 'hdi_rank_2021'])
energy.rename(columns={'Code': 'ISO'}, inplace=True)
energy.rename(columns={'Primary energy consumption per GDP (kWh/$)': 'Energy Intensity'}, inplace=True)
policies = policies[['country_iso', 'country', 'decision_date']]
hdi.columns = hdi.columns.str.replace('hdi_', '') #Remove 'hdi_' from column names

#Create row with country and year as index for datasets to be merged
gdp = gdp.melt(id_vars=['Country Name', 'ISO'], var_name='Year', value_name='GDP')
hdi = hdi.melt(id_vars=['country', 'ISO'], var_name='Year', value_name='HDI')

#Clean melted dataframes
gdp = gdp[gdp.Year != 'Unnamed: 65']
gdp = gdp[gdp.Year != 'Unnamed: 67']

# Merge dataframes with emissions
emissions['Year'] = emissions['Year'].astype(int)
gdp['Year'] = gdp['Year'].astype(int)
hdi['Year'] = hdi['Year'].astype(int)
energy['Year'] = energy['Year'].astype(int)
emissions_gdp = pd.merge(emissions, gdp, left_on=['ISO', 'Year'], right_on=['ISO', 'Year'])
emissions_gdp_hdi = pd.merge(emissions_gdp, hdi, left_on=['ISO', 'Year'], right_on=['ISO', 'Year'])
emissions_gdp_hdi_energy = pd.merge(emissions_gdp_hdi, energy, left_on=['ISO', 'Year'], right_on=['ISO', 'Year'], how='left')

# Remove duplicate columns
emissions_gdp_hdi_energy.drop(columns=['Country Name', 'country', 'Entity'], inplace=True)

# Remove sector breakdown columns and rename columns
emissions_gdp_hdi_energy.drop(columns=['Coal', 'Oil', 'Gas', 'Cement', 'Flaring', 'Other'], inplace=True)
emissions_gdp_hdi_energy = emissions_gdp_hdi_energy.rename(columns={'Total': 'Total CO2', 'Per Capita': 'CO2 Per Capita', 'GDP': 'GDP Per Capita'})

# Create Europeon Union dataset
eu = pd.DataFrame({'country': ['Austria', 'Belgium', 'Bulgaria', 'Croatia', 'Cyprus', 'Czech Republic', 'Denmark',
                              'Estonia', 'Finland', 'France', 'Germany', 'Greece', 'Hungary', 'Ireland', 'Italy',
                              'Latvia', 'Lithuania', 'Luxembourg', 'Malta', 'Netherlands', 'Poland', 'Portugal',
                               'Romania', 'Slovakia', 'Slovenia', 'Spain', 'Sweden'],
                   'ISO': ['AUT', 'BEL', 'BGR', 'HRV', 'CYP', 'CZE', 'DNK', 'EST', 'FIN', 'FRA', 'DEU', 'GRC', 'HUN',
                            'IRL', 'ITA', 'LVA', 'LTU', 'LUX', 'MLT', 'NLD', 'POL', 'PRT', 'ROU', 'SVK', 'SVN', 'ESP', 'SWE'], 
                   'Join Year': [1995, 1957, 2007, 2013, 2004, 2004, 1973, 2004, 1995, 1957, 1957, 1981, 2004, 1973, 1957, 2004,
                                  2004, 1957, 2004, 1958, 2004, 1986, 2007, 2004, 2004, 1986, 1995]})

# Create a dummy variable that indicates if a country is in the European Union for each year
emissions_gdp_hdi_energy = pd.merge(emissions_gdp_hdi_energy, eu, left_on='ISO', right_on='ISO', how='left')
emissions_gdp_hdi_energy['EU Member'] = np.where(emissions_gdp_hdi_energy['Year'] >= emissions_gdp_hdi_energy['Join Year'], 1, 0)
emissions_gdp_hdi_energy.drop(columns=['Join Year', 'country'], inplace=True)

# Create a secondary ISO code equal to 'EUE' in rows where EU Member is 1
emissions_gdp_hdi_energy['ISO_2'] = np.where(emissions_gdp_hdi_energy['EU Member'] == 1, 'EUE', np.nan)

# Merge policies dataset with emissions_gdp_hdi_energy dataset
emissions_gdp_hdi_energy_policy = pd.merge(emissions_gdp_hdi_energy, policies, left_on=['ISO', 'Year'],
                                           right_on=['country_iso', 'decision_date'], how='left')

# Create a dummy variable that indicates if a country had a climate policy ratified for each year
emissions_gdp_hdi_energy_policy['GHG Policy Ratified'] = np.where(emissions_gdp_hdi_energy_policy['country_iso'].notnull(), 1, 0)
emissions_gdp_hdi_energy_policy.drop(columns=['EU Member', 'country_iso', 'country', 'decision_date'], inplace=True)

# Repeat for EU countries
emissions_gdp_hdi_energy_policy = pd.merge(emissions_gdp_hdi_energy_policy, policies, left_on=['ISO_2', 'Year'],
                                     right_on=['country_iso', 'decision_date'], how='left')

# Set Policy to 1 for EU countries with decision_date
emissions_gdp_hdi_energy_policy['GHG Policy Ratified'] = np.where(emissions_gdp_hdi_energy_policy['country_iso'].notnull(), 1, emissions_gdp_hdi_energy_policy['GHG Policy Ratified'])
emissions_gdp_hdi_energy_policy.drop(columns=['country_iso', 'ISO_2', 'country', 'decision_date'], inplace=True)

# Remove duplicate rows
emissions_gdp_hdi_energy_policy = emissions_gdp_hdi_energy_policy.drop_duplicates()

# Set GRG policy ratified to 1 for any year after a ratified a policy 
emissions_gdp_hdi_energy_policy['GHG Policy Ratified'] = emissions_gdp_hdi_energy_policy.groupby('Country')['GHG Policy Ratified'].cummax()

#Create Temporary variables for HDI classifaction
emissions_gdp_hdi_energy_policy['Very High'] = 0
emissions_gdp_hdi_energy_policy['High'] = 0
emissions_gdp_hdi_energy_policy['Medium'] = 0
emissions_gdp_hdi_energy_policy['Low'] = 0

#Fill temporary variables with 1 if HDI is in the specified range
emissions_gdp_hdi_energy_policy['Very High'] = np.where(emissions_gdp_hdi_energy_policy['HDI'] > 0.8, 1, 0)
emissions_gdp_hdi_energy_policy['High'] = np.where((emissions_gdp_hdi_energy_policy['HDI'] <= 0.8) & (emissions_gdp_hdi_energy_policy['HDI'] > 0.7), 1, 0)  
emissions_gdp_hdi_energy_policy['Medium'] = np.where((emissions_gdp_hdi_energy_policy['HDI'] <= 0.7) & (emissions_gdp_hdi_energy_policy['HDI'] > 0.55), 1, 0)
emissions_gdp_hdi_energy_policy['Low'] = np.where(emissions_gdp_hdi_energy_policy['HDI'] <= 0.55, 1, 0)

#Create a new column that contains the HDI category for each country
emissions_gdp_hdi_energy_policy['HDI Category'] = 0
emissions_gdp_hdi_energy_policy['HDI Category'] = np.where(emissions_gdp_hdi_energy_policy['Very High'] == 1, 3, emissions_gdp_hdi_energy_policy['HDI Category'])
emissions_gdp_hdi_energy_policy['HDI Category'] = np.where(emissions_gdp_hdi_energy_policy['High'] == 1, 2, emissions_gdp_hdi_energy_policy['HDI Category'])
emissions_gdp_hdi_energy_policy['HDI Category'] = np.where(emissions_gdp_hdi_energy_policy['Medium'] == 1, 1, emissions_gdp_hdi_energy_policy['HDI Category'])
emissions_gdp_hdi_energy_policy['HDI Category'] = np.where(emissions_gdp_hdi_energy_policy['Low'] == 1, 0, emissions_gdp_hdi_energy_policy['HDI Category'])

#Remove temporary variables
emissions_gdp_hdi_energy_policy = emissions_gdp_hdi_energy_policy.drop(columns=['Very High', 'High', 'Medium', 'Low'])

# Save cleaned data
emissions_gdp_hdi_energy_policy = emissions_gdp_hdi_energy_policy[['Country', 'ISO', 'Year', 'Total CO2', 'CO2 Per Capita', 'GDP Per Capita', 'HDI', 'HDI Category', 'Energy Intensity', 'GHG Policy Ratified']]

emissions_gdp_hdi_energy_policy.head()

# Function that creates a list of all countries with missing data in a specified column
def missing_data(df, column) -> list:
    """Return a list of countries with missing data in the specified column"""
    missing = df[df[column].isnull()]['Country'].unique()
    return missing

# Create subset of data with no observations past 2017
e2017 = emissions_gdp_hdi_energy_policy[emissions_gdp_hdi_energy_policy['Year'] <= 2017]

# Number of Obersevations before data cleaning
before = emissions_gdp_hdi_energy_policy['GDP Per Capita'].count()

# Create lists of countries with missing data
m_gdp = missing_data(emissions_gdp_hdi_energy_policy, 'GDP Per Capita')
m_hdi = missing_data(emissions_gdp_hdi_energy_policy, 'HDI')
m_energy = missing_data(e2017, 'Energy Intensity')

# Merge lists of countries with missing data into one list with no duplicates
missing = list(set(m_gdp) | set(m_hdi) | set(m_energy))

# Fuction that removes all observations with specified country name
def remove_country(df, country) -> pd.DataFrame:
    """Return a dataframe with all observations for the specified country removed"""
    df = df[df['Country'] != country]
    return df   

# Remove countries with missing data
for country in missing:
    emissions_gdp_hdi_energy_policy = remove_country(emissions_gdp_hdi_energy_policy, country)

# Number of Obersevations after data cleaning
after = emissions_gdp_hdi_energy_policy['GDP Per Capita'].count()

# Create a new variable for clean called ID that is the ISO code and year
emissions_gdp_hdi_energy_policy['ID'] = emissions_gdp_hdi_energy_policy['ISO'] + emissions_gdp_hdi_energy_policy['Year'].astype(str)

# Save cleaned data
emissions_gdp_hdi_energy_policy.to_csv('cleaned_data.csv', index=False)

clean = pd.read_csv('cleaned_data.csv')

# Print summary of data cleaning
print(f'Number of observations dropped from data cleaning: {before - after}')
print(f'Number of countries dropped from data cleaning: {len(missing)}')
print(f'Number of observations in cleaned data: {after}')
clean.head()

Number of observations dropped from data cleaning: 2263
Number of countries dropped from data cleaning: 84
Number of observations in cleaned data: 3520

clean.drop(columns=['Year']).describe()

# Create time series plot that show the total number of countries with a GHG policy ratified from 1990-2021
policy_time_series = clean.groupby('Year')['GHG Policy Ratified'].sum().reset_index()

# Create bar graph that plots total CO2 emissions produced each year by countries with a GHG policy ratified
# Group all countries by year and sum the total CO2 emissions for each year
total_emissions = clean.groupby('Year')['Total CO2'].sum().reset_index()
# Sum the total CO2 emissions for each year only including countries with a GHG policy ratified
policy_emissions = clean[clean['GHG Policy Ratified'] == 1].groupby('Year')['Total CO2'].sum().reset_index()
# Merge the two dataframes on year
emissions_policy = pd.merge(total_emissions, policy_emissions, on='Year')

# Create a figure with two subplots
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(8, 8))

# Plot the first graph in the first subplot
ax1.bar(policy_time_series['Year'], policy_time_series['GHG Policy Ratified'], color='seagreen')
ax1.set_title('Total Number Countries with any GHG Policy Ratified')
ax1.set_xlabel('Year')
ax1.set_ylabel('Total Countries')
ax1.spines['top'].set_visible(False)
ax1.spines['right'].set_visible(False)
ax1.yaxis.grid(color='lightgrey')


# Plot the second graph in the second subplot
ax2.bar(emissions_policy['Year'], emissions_policy['Total CO2_x'], label='Total CO2 Emissions', color='tomato')
ax2.bar(emissions_policy['Year'], emissions_policy['Total CO2_y'], label='Emissions Produced by Countries with GHG targets')
ax2.set_title('Total CO2 Emissions Produced Each Year')
ax2.set_xlabel('Year')
ax2.set_ylabel('Total CO2 Emissions (MtCO2)')
ax2.spines['top'].set_visible(False)
ax2.spines['right'].set_visible(False)
ax2.yaxis.grid(color='lightgrey')
ax2.legend()

# Adjust the spacing between subplots
fig.tight_layout()

# Show the figure
plt.show()

# Create new variable Carbon Intensity Per Capita
clean['Carbon Intensity Per Capita'] = clean['CO2 Per Capita'] / clean['Energy Intensity']

# Create Copy of cleaned data
clean_pc_5 = clean.copy()

# Function that calculates the percent change over the past X years in a specified column for each country and appends the results to the dataframe
def percent_change(df: pd.DataFrame, group: str, column: str, time:int=5) -> pd.DataFrame:
    """Calculates the percent change over the past X years in a specified column for each country and appends the results to the dataframe"""
    df[f'{column} Percent Change'] = df.groupby(group)[column].pct_change(periods=time)
    return df

emissions_pc = percent_change(clean_pc_5, 'Country', 'CO2 Per Capita')
emissions_pc = percent_change(emissions_pc, 'Country', 'GDP Per Capita')
emissions_pc = percent_change(emissions_pc, 'Country', 'Energy Intensity')
emissions_pc = percent_change(emissions_pc, 'Country', 'Carbon Intensity Per Capita')

# Totals of policy ratified and non-ratified countries in each quadrant in a given year
def four_quad_totals(df: pd.DataFrame, column_1: str, column_2:str, binary_c: str, year:int) -> tuple[int]:
    """Take a dataframe and two column names with float values. Returns a list of tuples that contain the count of
    countries with postive values in both columns, negative values in both columns, positive in column_1 negative in column_2, and
    nagative in column_1 positive in column_2 for a specified year. The first index of the tuple is the count of countires with 
    binary_c equal to 0 and the second index of the tuple is the count of countries with binary_c equal to 1"""
    # Both columns are positive and binary_c is 0
    pos_pos_0 = df[(df[column_1] > 0) & (df[column_2] > 0) & (df[binary_c] == 0) & (df['Year'] == year)][column_1].count()
    # Both columns are positive and binary_c is 1
    pos_pos_1 = df[(df[column_1] > 0) & (df[column_2] > 0) & (df[binary_c] == 1) & (df['Year'] == year)][column_1].count()
    # Both columns are negative and binary_c is 0
    neg_neg_0 = df[(df[column_1] < 0) & (df[column_2] < 0) & (df[binary_c] == 0) & (df['Year'] == year)][column_1].count()
    # Both columns are negative and binary_c is 1
    neg_neg_1 = df[(df[column_1] < 0) & (df[column_2] < 0) & (df[binary_c] == 1) & (df['Year'] == year)][column_1].count()
    # column_1 is positive and column_2 is negative and binary_c is 0
    pos_neg_0 = df[(df[column_1] > 0) & (df[column_2] < 0) & (df[binary_c] == 0) & (df['Year'] == year)][column_1].count()
    # column_1 is positive and column_2 is negative and binary_c is 1
    pos_neg_1 = df[(df[column_1] > 0) & (df[column_2] < 0) & (df[binary_c] == 1) & (df['Year'] == year)][column_1].count()
    # column_1 is negative and column_2 is positive and binary_c is 0
    neg_pos_0 = df[(df[column_1] < 0) & (df[column_2] > 0) & (df[binary_c] == 0) & (df['Year'] == year)][column_1].count()
    # column_1 is negative and column_2 is positive and binary_c is 1
    neg_pos_1 = df[(df[column_1] < 0) & (df[column_2] > 0) & (df[binary_c] == 1) & (df['Year'] == year)][column_1].count()
    return (pos_pos_0, pos_pos_1), (neg_neg_0, neg_neg_1), (pos_neg_0, pos_neg_1), (neg_pos_0, neg_pos_1)


# Function that creates multiple scatter plots that show the relationship between two columns for every country in specified years
def multiple_scatter_plots(df: pd.DataFrame, time: str, years: list[int], x: str, y: str, xrange: list[int], yrange: list[int], log_x: bool = False, log_y: bool = False) -> None:
    """Creates multiple scatter plots that show the relationship between two columns for every country in specified years"""
    fig, axs = plt.subplots(1, len(years), figsize=(15, 5), sharey=True)  # Add sharey=True to share the y-axis
    for i, year in enumerate(years):
        year_data = df[df['Year'] == year]
        if log_x:
            axs[i].set_xscale('log')
        if log_y:
            axs[i].set_yscale('log')
        axs[i].scatter(year_data[x], year_data[y], c=year_data['GHG Policy Ratified'].apply(lambda x: 'green' if x == 1 else 'red'), s=7)
        axs[i].set_xlabel(x)
        if year == 4:
            axs[i].set_title(f'Change over 4 years prior to GHG policy ratified')
        elif year == 8:
            axs[i].set_title(f'Change over 4 years after GHG policy ratified')
        else:
            axs[i].set_title(f'Change over last {time} years, {year}')
        
        axs[i].spines['top'].set_visible(False)
        axs[i].spines['right'].set_visible(False)
        axs[i].axvline(x=0, linestyle='dotted', color='gray')
        axs[i].axhline(y=0, linestyle='dotted', color='gray')
        axs[i].set_xlim(left=xrange[0], right=xrange[1])
        axs[i].set_ylim(bottom=yrange[0], top=yrange[1])
        both_pos, both_neg, gdppos_co2neg, gdpneg_co2pos = four_quad_totals(df, x, y, 'GHG Policy Ratified', year)
        # positions
        btm_l = [both_neg, 0.05, 0.05, 'left', 'bottom']
        btm_r = [gdppos_co2neg, 0.95, 0.05, 'right', 'bottom']
        top_l = [gdpneg_co2pos, 0.05, 0.95, 'left', 'top']
        top_r = [both_pos, 0.95, 0.95, 'right', 'top']
        corners = [btm_l, btm_r, top_l, top_r]
        # Textboxs
        for j in corners:
            axs[i].text(j[1], j[2], f"GHG Policy:\n{j[0][1]}/{j[0][1] + j[0][0]}", horizontalalignment=j[3],
                         verticalalignment=j[4], transform=axs[i].transAxes)
    axs[0].set_ylabel(y)
    plt.show()

# report_outliers(emissions_pc, 'GDP Per Capita Percent Change', 'CO2 Per Capita Percent Change', year, [-0.5, 1], [-0.5, 1.25])
multiple_scatter_plots(emissions_pc, '5', [1996, 2001, 2006], 'GDP Per Capita Percent Change', 'Carbon Intensity Per Capita Percent Change', [-1, 1], [-1, 1])
multiple_scatter_plots(emissions_pc, '5', [2011, 2016, 2021], 'GDP Per Capita Percent Change', 'Carbon Intensity Per Capita Percent Change', [-1, 1], [-1, 1])

# Create a dataset with only observations with a GHG policy ratified
policy = clean[clean['GHG Policy Ratified'] == 1]
# Create a dataset with only observations with no GHG policy ratified
no_policy = clean[clean['GHG Policy Ratified'] == 0]

# Create dataframes that contain the first observation for each country in policy
first_policy = policy.groupby('Country').first().reset_index()
first_policy = first_policy[['ISO', 'Year']]
first_policy.rename(columns={'Year': 'GHG Policy Year'}, inplace=True)

# Merge first_policy with policy and no_policy
policy = pd.merge(policy, first_policy, on='ISO', how='left')
no_policy = pd.merge(no_policy, first_policy, on='ISO', how='left')

val = 4

# Remove all observations that are not within the window
policy = policy[policy['Year'] - policy['GHG Policy Year'] <= val-1]
no_policy = no_policy[no_policy['GHG Policy Year'] - no_policy['Year'] <= val]

# Remove countries with missing observations
policy = policy.groupby('Country').filter(lambda x: len(x) == val)
no_policy = no_policy.groupby('Country').filter(lambda x: len(x) == val)

# Remove all countries with missing data in window
missing_policy_ci = missing_data(policy, 'Carbon Intensity Per Capita')
missing_policy_gdp = missing_data(policy, 'GDP Per Capita')
missing_no_policy_ci = missing_data(no_policy, 'Carbon Intensity Per Capita')
missing_no_policy_gdp = missing_data(no_policy, 'GDP Per Capita')
missing = list(set(missing_policy_ci) | set(missing_policy_gdp) | set(missing_no_policy_ci) | set(missing_no_policy_gdp))
for country in missing:
    policy = remove_country(policy, country)
    no_policy = remove_country(no_policy, country)

# Create dataset that removes all observations from clean that are not in policy or no_policy based on ID
clean_no_policy = clean[clean['ID'].isin(policy['ID']) | clean['ID'].isin(no_policy['ID'])]

# Remove all countries with less than 8 observations
clean_no_policy = clean_no_policy.groupby('Country').filter(lambda x: len(x) == 8)

# Duplicate 'Year' column
clean_no_policy['Original Year'] = clean_no_policy['Year']

# Change values in 'Year' column to repeat 1-8
clean_no_policy['Year'] = clean_no_policy.groupby('Country').cumcount() + 1

emissions_pc_4 = percent_change(clean_no_policy, 'Country', 'GDP Per Capita', 3)
emissions_pc_4 = percent_change(emissions_pc_4, 'Country', 'Carbon Intensity Per Capita', 3)

multiple_scatter_plots(emissions_pc_4, 4, [4, 8], 'GDP Per Capita Percent Change', 'Carbon Intensity Per Capita Percent Change', [-1, 1], [-1, 1])

# Create dataframe with only 1990 observations
clean_1990 = clean[clean['Year'] == 1990]
# Create dataframe with only 2018 observations
clean_2018 = clean[clean['Year'] == 2018]
# Join the two dataframes on country
change = pd.merge(clean_1990, clean_2018, on='Country', suffixes=('_1990', '_2018'))

# Create variable that shows the change in Carbon intensity per capita, CO2 per capita, HDI, GDP Per capita, and Energy Intensity
change['Carbon Intensity Change'] = change['Carbon Intensity Per Capita_2018'] - change['Carbon Intensity Per Capita_1990']
change['HDI Change'] = change['HDI_2018'] - change['HDI_1990']
change['CO2 Change'] = change['CO2 Per Capita_2018'] - change['CO2 Per Capita_1990']
change['GDP Change'] = change['GDP Per Capita_2018'] - change['GDP Per Capita_1990']
change['Energy Intensity Change'] = change['Energy Intensity_2018'] - change['Energy Intensity_1990']

# Calculate percent change in Carbon intensity per capita, CO2 per capita, HDI, GDP Per capita, and Energy Intensity
change['Carbon Intensity PC'] = change['Carbon Intensity Change'] / change['Carbon Intensity Per Capita_1990']
change['HDI PC'] = change['HDI Change'] / change['HDI_1990']
change['CO2 PC'] = change['CO2 Change'] / change['CO2 Per Capita_1990']
change['GDP PC'] = change['GDP Change'] / change['GDP Per Capita_1990']
change['Energy Intensity PC'] = change['Energy Intensity Change'] / change['Energy Intensity_1990']

# Create variable ln GDP change that is the natural log of the GDP change
change['ln GDP Change'] = np.log(change['GDP Change'] + 1)


# make a copy of change dataframe called change_1
change_1 = change.copy()

# Remove all contrieis with High or Very High HDI in 1990
change_h = change[change['HDI_1990'] < 0.7]
# Remove all countries with Low or Medium HDI in 2018
change_l = change_h[change_h['HDI_2018'] > 0.7]
plt.figure(figsize=(10, 6))
plt.scatter(change['GDP Change'], change['CO2 Change'], color='green', s=6)
plt.scatter(change_h['GDP Change'], change_h['CO2 Change'], color='blue', s=6)
plt.scatter(change_l['GDP Change'], change_l['CO2 Change'], color='red', s=6)
plt.gca().spines['top'].set_visible(False)
plt.gca().spines['right'].set_visible(False)
plt.gca().set_ylim(bottom=-10, top=10)
plt.gca().set_xlim(left=0, right=45000)
plt.gca().axhline(y=0, linestyle='dotted', color='gray')
plt.title('Change in GDP per Capita vs. Change in CO2 per Capita, 1990-2018')
plt.xlabel('Change in GDP per Capita ($)')
plt.ylabel('Change in CO2 per Capita (MtCO2)')
plt.text(100, -4, 'Countries with HDI less than .7', color='blue')
plt.text(2500, -4.7, 'from 1990-2018', color='blue')
plt.plot([6000, 2000], [-3.5, -.1], color='blue')
plt.text(10000, 8.5, 'Countries with HDI less than .7 in 1990', color='red')
plt.text(11000, 7.8, 'and HDI greater than .7 in 2018', color='red')
plt.plot([10300, 16000], [1.25, 7.3], color='red')
plt.plot([18500, 16000], [2.6, 7.3], color='red')
plt.text(27000, -7, 'Countries with HDI greater than .7', color='green')
plt.text(30800, -7.6, 'from 1990-2018', color='green')
plt.plot([34800, 34000], [-6.5, -4.5], color='green')
plt.gca().yaxis.grid(color='lightgrey')
plt.show()

# Create a geopandas dataframe with the world map
world = gpd.read_file(gpd.datasets.get_path("naturalearth_lowres"))
world = world.set_index("iso_a3")
# Create dataframe with only observations for 1990 and 2018
clean_1990 = clean[clean['Year'] == 1990]
clean_2018 = clean[clean['Year'] == 2018]
warnings.filterwarnings("ignore")
# Create a geopandas dataframe with the world map
world_2018 = world.join(clean_2018.set_index('ISO'), how='outer', rsuffix='_2018')
world_1990 = world.join(clean_1990.set_index('ISO'), how='outer', rsuffix='_1990')

def plot_2_maps(geo_df_left, geo_df_right, column, min_val, max_val, unit=0, color='Reds'):
    """Plot two maps side by side for the specified column"""
    fig, axs = plt.subplots(1, 2, figsize=(15, 7))
    fig.suptitle(column + ' by Country, 1990 vs 2018', fontsize=20)
    fig.set_facecolor((0.96, 0.96, 0.96))
    # Create a single legend for both subplots
    cax = fig.add_axes([0.49, 0.2, 0.62, 0.2])
    # Plot 1990 map
    geo_df_left.plot(column=column, cmap=color, linewidth=0.8, ax=axs[0], edgecolor='0.8', legend=False, vmin=min_val, vmax=max_val, missing_kwds={"color": "lightgrey", "label": "Missing values"})
    axs[0].set_title(column + ' in 1990', fontsize=15)
    axs[0].spines['top'].set_visible(False)
    axs[0].spines['right'].set_visible(False)
    axs[0].set_xlabel('Longitude')
    axs[0].set_ylabel('Latitude')
    # Plot 2018
    geo_df_right.plot(column=column, cmap=color, linewidth=0.8, ax=axs[1], edgecolor='0.8', legend=False, vmin=min_val, vmax=max_val, missing_kwds={"color": "lightgrey", "label": "Missing values"})
    axs[1].set_title(column + ' in 2018', fontsize=15)
    axs[1].spines['top'].set_visible(False)
    axs[1].spines['right'].set_visible(False)
    axs[1].set_xlabel('Longitude')
    # Create a single legend for both subplots
    sm = plt.cm.ScalarMappable(cmap=color, norm=plt.Normalize(vmin=min_val, vmax=max_val))
    sm.set_array([])
    fig.colorbar(sm, cax=cax)
    if unit == 0:
        fig.text(0.5, 0.84, column, ha='center')
    else:
        fig.text(0.5, 0.87, column, ha='center')
        fig.text(0.5, 0.84,' (' + unit + ')', ha='center')
    fig.sharey=True
    plt.show()

def plot_2_mapsh(geo_df_left, geo_df_right, column, min_val, max_val, unit=0, color='Reds', start='1990', end='2018'):
    """Plot two maps stacked on top of each other for the specified column"""
    fig, axs = plt.subplots(2, 1, figsize=(25, 20))
    fig.suptitle(column + ' by Country, ' + start + ' vs ' + end, fontsize=20)
    fig.set_facecolor((0.96, 0.96, 0.96))
    # Create a single legend for both subplots
    cax = fig.add_axes([.36, 0.485, 0.28, .02])
    # Plot 1990 map
    geo_df_left.plot(column=column, cmap=color, linewidth=0.8, ax=axs[0], edgecolor='0.8', legend=False, vmin=min_val, vmax=max_val, missing_kwds={"color": "lightgrey", "label": "Missing values"})
    axs[0].set_title(column + ' in ' + start, fontsize=15)
    axs[0].spines['top'].set_visible(False)
    axs[0].spines['right'].set_visible(False)
    axs[0].set_xlabel('Longitude')
    axs[0].set_ylabel('Latitude')
    # Plot 2018
    geo_df_right.plot(column=column, cmap=color, linewidth=0.8, ax=axs[1], edgecolor='0.8', legend=False, vmin=min_val, vmax=max_val, missing_kwds={"color": "lightgrey", "label": "Missing values"})
    axs[1].set_title(column + ' in ' + end, fontsize=15)
    axs[1].spines['top'].set_visible(False)
    axs[1].spines['right'].set_visible(False)
    axs[1].set_xlabel('Longitude')
    # Create a single legend for both subplots
    sm = plt.cm.ScalarMappable(cmap=color, norm=plt.Normalize(vmin=min_val, vmax=max_val))
    sm.set_array([])
    fig.colorbar(sm, cax=cax, orientation="horizontal")
    if unit == 0:
        fig.text(0.675, 0.49, column, ha='center')
    else:
        fig.text(0.675, 0.495, column, ha='center')
        fig.text(0.675, 0.485,' (' + unit + ')', ha='center')
    fig.sharex=True
    fig.tight_layout    
    plt.show()

plot_2_mapsh(world_1990, world_2018, 'CO2 Per Capita', 0, 25, 'MtCO2')
plot_2_mapsh(world_1990, world_2018, 'GDP Per Capita', 0, 80000, '$', 'Blues')

plot_2_mapsh(world_1990, world_2018, 'HDI', .2, 1, 0, 'RdYlGn')
plot_2_mapsh(world_1990, world_2018, 'Carbon Intensity Per Capita', 0, 10, 'MtCO2/kWh/$', 'RdYlGn_r')

web_url = 'https://www.carbonbrief.org/analysis-which-countries-are-historically-responsible-for-climate-change/'
headers = {'User-Agent': 'Mozilla/5.0 (Windows NT 10.0; Win64; x64) AppleWebKit/537.36 (KHTML, like Gecko) Chrome/58.0.3029.110 Safari/537.3'}
response = requests.get(web_url, headers=headers)

soup_object = BeautifulSoup(response.content)
data_table = soup_object.find_all('table')[0]

ix = 0
cumulative_emissions = pd.DataFrame(columns = ['rank', 'country', 'cumulative emissions per population 2021',
                                                'rank_0', 'country_0', 'cumulative emissions per capita 2021'])

all_values = data_table.find_all('tr')
for row in all_values[1:]:
    values = row.find_all('td')
    rank = values[0].text.strip()
    country = values[1].text.strip()
    cumulative_per_population = values[2].text.strip()
    rank_0 = values[3].text.strip()
    country_0 = values[4].text.strip()
    cumulative_per_capita = values[5].text.strip()
    cumulative_emissions.loc[ix] = [rank, country, cumulative_per_population, rank_0, country_0, cumulative_per_capita]
    ix += 1

cumulative_emissions.head()

# Split Cumlative emissions into two dataframes
cumulative_per_population = cumulative_emissions[['country', 'rank', 'cumulative emissions per population 2021']]
cumulative_per_capita = cumulative_emissions[['country_0', 'rank_0', 'cumulative emissions per capita 2021']]
cumulative_per_population.rename(columns={'country': 'Country', 'rank': 'Rank Per Population', 'cumulative emissions per population 2021': 'Cumulative Emissions Per Population 2021'}, inplace=True)
cumulative_per_capita.rename(columns={'country_0': 'Country', 'rank_0': 'Rank Per Capita', 'cumulative emissions per capita 2021': 'Cumulative Emissions Per Capita 2021'}, inplace=True)

# Create data with only 2021 observations from clean
clean_2021 = clean[clean['Year'] == 2021]

# Merge clean with cumulative_per_population
population_merged = pd.merge(cumulative_per_population, clean_2021, on='Country', how='inner')
per_capita_merged = pd.merge(cumulative_per_capita, clean_2021, on='Country', how='inner')

# Create a scatter plot for 2 columns 
plt.figure(figsize=(10, 6))
plt.scatter(population_merged['HDI'], population_merged['Cumulative Emissions Per Population 2021'], color='green')
plt.gca().spines['top'].set_visible(False)
plt.gca().spines['right'].set_visible(False)
plt.gca().invert_yaxis()
# plt.gca().set_ylim(bottom=-10, top=10)
plt.gca().set_xlim(left=.55, right=1)
plt.title('Cumulative Emissions Per Population 2021 vs. HDI 2021')
plt.xlabel('HDI in 2021')
plt.ylabel('Cumulative Emissions Per Population 2021')
plt.axvline(x=.7, color='grey', linestyle='--')
plt.show()
plt.figure(figsize=(10, 6))
plt.scatter(per_capita_merged['HDI'], per_capita_merged['Cumulative Emissions Per Capita 2021'], color='green')
plt.gca().spines['top'].set_visible(False)
plt.gca().spines['right'].set_visible(False)
plt.gca().invert_yaxis()
# plt.gca().set_ylim(bottom=-10, top=10)
plt.gca().set_xlim(left=.55, right=1)
plt.title('Cumulative Emissions Per Capita 2021 vs. HDI 2021')
plt.xlabel('HDI in 2021')
plt.ylabel('Cumulative Emissions Per Capita 2021')
plt.axvline(x=.7, color='grey', linestyle='--')
plt.show()

# Put temp data into a dataframe
temp = pd.read_csv('GlobalLandTemperaturesByCountry.csv')
temp['dt'] = pd.to_datetime(temp['dt'])
# Remove all observation from before 1990
temp = temp[temp['dt'].dt.year >= 1990]

# Create a new column that contains the average temperature for each country for each year
temp['Year'] = temp['dt'].dt.year
temp = temp.groupby(['Country', 'Year'])['AverageTemperature'].mean().reset_index()

# Rename countries to match ISO codes
temp['Country'] = temp['Country'].replace({'Congo (Democratic Republic Of The)': 'Congo (Kinshasa)', 'Congo': 'Congo (Brazzaville)', 'Côte D\'Ivoire': 'Ivory Coast', 'Myanmar': 'Burma', 'Réunion': 'Reunion', 'São Tomé And Príncipe': 'Sao Tome and Principe', 'South Georgia And The South Sandwich Isla': 'South Georgia and the South Sandwich Islands', 'St. Barthélemy': 'Saint Barthelemy', 'St. Vincent And The Grenadines': 'Saint Vincent and the Grenadines', 'Syrian Arab Republic': 'Syria', 'Taiwan': 'Taiwan*', 'Tanzania': 'Tanzania*', 'United States': 'USA', 'Venezuela': 'Venezuela*', 'Vietnam': 'Vietnam*'})

# Merge temperature data with clean data
clean_temp = pd.merge(clean, temp, on=['Country', 'Year'], how='left')

clean_temp

# Create dataframes that only contain data for 1990 and 2013
temp_1990 = clean_temp[clean_temp['Year'] == 1990]
temp_2013 = clean_temp[clean_temp['Year'] == 2013]


# Create a geopandas dataframe with the world map
world_temp = world.join(temp_1990.set_index('ISO'), how='outer', rsuffix='_1990')
world_temp = world_temp.join(temp_2013.set_index('ISO'), how='outer', rsuffix='_2013')

# Find the percent increase in temperature from 1990 to 2013
world_temp['Temperature Change'] = world_temp['AverageTemperature_2013'] - world_temp['AverageTemperature']
world_temp['Temperature Percent Change'] = world_temp['Temperature Change'] / world_temp['AverageTemperature']

# Find the percent increase in GDP per Capita from 1960 to 2013
world_temp['GDP Change'] = world_temp['GDP Per Capita_2013'] - world_temp['GDP Per Capita']
world_temp['GDP Percent Change'] = world_temp['GDP Change'] / world_temp['GDP Per Capita']

# Find percent increase in CO2 per Capita from 1960 to 2013
world_temp['CO2 Change'] = world_temp['CO2 Per Capita_2013'] - world_temp['CO2 Per Capita']
world_temp['CO2 Percent Change'] = world_temp['CO2 Change'] / world_temp['CO2 Per Capita']

# Find percent increase in Energy Intensity from 1960 to 2013
world_temp['Energy Intensity Change'] = world_temp['Energy Intensity_2013'] - world_temp['Energy Intensity']
world_temp['Energy Intensity Percent Change'] = world_temp['Energy Intensity Change'] / world_temp['Energy Intensity']

# The change in HDI from 1990 to 2013
world_temp['HDI Change'] = world_temp['HDI_2013'] - world_temp['HDI']
world_temp['HDI Percent Change'] = world_temp['HDI Change'] / world_temp['HDI']


def plot_map(geo_df_left, column, min_val, max_val, year, title, unit=0, color='Reds'):
    """Plot a map for the specified column"""
    fig, ax = plt.subplots(1, 1, figsize=(12, 6))
    fig.suptitle(title, fontsize=20)
    fig.set_facecolor((0.96, 0.96, 0.96))
    # Create a single legend for both subplots
    cax = fig.add_axes([0.93, 0.2, 0.05, 0.63])
    # Plot map
    geo_df_left.plot(column=column, cmap=color, linewidth=0.8, ax=ax, edgecolor='0.8', legend=False, vmin=min_val, vmax=max_val, missing_kwds={"color": "lightgrey", "label": "Missing values"})
    ax.spines['top'].set_visible(False)
    ax.spines['right'].set_visible(False)
    ax.set_xlabel('Longitude')
    ax.set_ylabel('Latitude')
    # Create a single legend for both subplots
    sm = plt.cm.ScalarMappable(cmap=color, norm=plt.Normalize(vmin=min_val, vmax=max_val))
    sm.set_array([])
    fig.colorbar(sm, cax=cax)
    if unit == 0:
        fig.text(0.95, 0.85, column, ha='center')
    else:
        fig.text(0.95, 0.88, column, ha='center')
        fig.text(0.95, 0.85,' (' + unit + ')', ha='center')
    plt.show()

plot_map(world_temp, 'Temperature Change', -4, 4, 2013, 'Change in Average Temperature by country from 1990 to 2013', 0, 'RdYlBu_r')

clean_temp = percent_change(clean_temp, 'Country', 'AverageTemperature', 3)
clean_temp = percent_change(clean_temp, 'Country', 'Carbon Intensity Per Capita', 3)

multiple_scatter_plots(clean_temp, '3', [1994, 1997, 2000], 'Carbon Intensity Per Capita Percent Change', 'AverageTemperature Percent Change', [-.5, 1], [-.35, .35])
multiple_scatter_plots(clean_temp, '3', [2003, 2006, 2009], 'Carbon Intensity Per Capita Percent Change', 'AverageTemperature Percent Change', [-.5, 1], [-.35, .35])

# Regression libraries
from stargazer.stargazer import Stargazer
from IPython.core.display import HTML
import statsmodels.api as sm
import statsmodels.formula.api as smf
# .png display package
from IPython.display import Image

# Open cleaned_data.csv
reg_df = pd.read_csv('cleaned_data.csv')

# Rename columns removing capitals and spaces
reg_df.rename(columns = {'Country':'country', 'Year':'year', 'CO2 Per Capita':'co2_pc', 'GDP Per Capita':'gdp_pc', 'HDI':'hdi', 'Energy Intensity':'ei', 'GHG Policy Ratified':'ghg', 'HDI Category':'hdi_cat'}, inplace = True)

# # Create Carbon Intensity Per Capita
reg_df['ci_pc'] = reg_df['co2_pc'] / reg_df['ei']

# Create a new column that number each country from 0 to 109 (110 countries)
reg_df['country_fe'] = reg_df['country'].astype('category')
reg_df['country_fe'] = reg_df['country_fe'].cat.codes

# Create variable that has ln(GDP per Capita)
reg_df['ln_gdp_pc'] = np.log(reg_df['gdp_pc'])

# Create 'Change in CO2 Per Capita' and 'Change in GDP Per Capita' columns
reg_df['change_co2'] = reg_df.groupby('country')['co2_pc'].diff()
reg_df['change_gdp'] = reg_df.groupby('country')['gdp_pc'].diff()
reg_df['change_ln_gdp'] = reg_df.groupby('country')['ln_gdp_pc'].diff()
reg_df['change_hdi'] = reg_df.groupby('country')['hdi'].diff()
reg_df['change_ei'] = reg_df.groupby('country')['ei'].diff()

# Drop 1990 observations
reg_df = reg_df[reg_df['year'] != 1990]

# Make 'Country FE' a categorical variable
reg_df['country_FE'] = reg_df['country_fe'].astype('category')

# Create clean_e with observations prior to 2019 for energy intensity regressions
reg_df = reg_df[reg_df['year'] <= 2018]

warnings.filterwarnings("ignore")

# Regression models, loggdp
reg = smf.ols(formula='change_co2 ~ change_ln_gdp + C(country_FE)', data=reg_df)
reg2 = smf.ols(formula='change_co2 ~ change_ln_gdp + ghg + C(country_FE)', data=reg_df)
reg3 = smf.ols(formula='change_co2 ~ change_ln_gdp + ghg + change_ei + C(country_FE)', data=reg_df)
reg4 = smf.ols(formula='change_co2 ~ change_ln_gdp + ghg + change_ei + C(hdi_cat) + C(country_FE)', data=reg_df)

# Fit the regression models
ares = reg.fit(cov_type='cluster', cov_kwds={'groups': reg_df['country_fe']})
ares2 = reg2.fit(cov_type='cluster', cov_kwds={'groups': reg_df['country_fe']})
ares3 = reg3.fit(cov_type='cluster', cov_kwds={'groups': reg_df['country_fe']})
ares4 = reg4.fit(cov_type='cluster', cov_kwds={'groups': reg_df['country_fe']})

# Create Stargazer object
stargazer = Stargazer([ares, ares2, ares3, ares4])

# print(stargazer.render_latex())

# Load image of regression results table without Country Fixed Effects coefficients 
Image(filename='reg_table1.png')

def set_base_year(df, year):
    """Return the dataframe with a new variable called 'HDI Category ' + year that has a country's HDI 
    category for the specified year for all observations."""
    temp_hdi = df[df['year'] == year][['country', 'hdi_cat']]
    temp_hdi.rename(columns={'hdi_cat': 'hdi_cat_' + str(year)}, inplace=True)
    base = pd.merge(df, temp_hdi, on='country', how='left')
    return base

base2015 = set_base_year(reg_df, 2015)

# Create 4 new dataframes, one for countries with Low HDI, one for countries with Medium HDI, one for countries with High HDI, and one for countries with Very High HDI
low_hdi = base2015[base2015['hdi_cat_2015'] == 0]
medium_hdi = base2015[base2015['hdi_cat_2015'] == 1]
high_hdi = base2015[base2015['hdi_cat_2015'] == 2]
very_high_hdi = base2015[base2015['hdi_cat_2015'] == 3]

# Regression models, low_hdi
reg = smf.ols(formula='change_co2 ~ change_ln_gdp + C(country_FE)', data=low_hdi)
reg2 = smf.ols(formula='change_co2 ~ change_ln_gdp + ghg + C(country_FE)', data=low_hdi)
reg3 = smf.ols(formula='change_co2 ~ change_ln_gdp + ghg + change_ei + C(country_FE)', data=low_hdi)

# Fit the regression models
res = reg.fit(cov_type='cluster', cov_kwds={'groups': low_hdi['country_fe']})
res2 = reg2.fit(cov_type='cluster', cov_kwds={'groups': low_hdi['country_fe']})
res3 = reg3.fit(cov_type='cluster', cov_kwds={'groups': low_hdi['country_fe']})

# Regression models, medium_hdi
reg4 = smf.ols(formula='change_co2 ~ change_ln_gdp + C(country_FE)', data=medium_hdi)
reg5 = smf.ols(formula='change_co2 ~ change_ln_gdp + ghg + C(country_FE)', data=medium_hdi)
reg6 = smf.ols(formula='change_co2 ~ change_ln_gdp + ghg + change_ei + C(country_FE)', data=medium_hdi)

# Fit the regression models
res4 = reg4.fit(cov_type='cluster', cov_kwds={'groups': medium_hdi['country_fe']})
res5 = reg5.fit(cov_type='cluster', cov_kwds={'groups': medium_hdi['country_fe']})
res6 = reg6.fit(cov_type='cluster', cov_kwds={'groups': medium_hdi['country_fe']})

# Regression models, high_hdi
reg7 = smf.ols(formula='change_co2 ~ change_ln_gdp + C(country_FE)', data=high_hdi)
reg8 = smf.ols(formula='change_co2 ~ change_ln_gdp + ghg + C(country_FE)', data=high_hdi)
reg9 = smf.ols(formula='change_co2 ~ change_ln_gdp + ghg + change_ei + C(country_FE)', data=high_hdi)

# Fit the regression models
res7 = reg7.fit(cov_type='cluster', cov_kwds={'groups': high_hdi['country_fe']})
res8 = reg8.fit(cov_type='cluster', cov_kwds={'groups': high_hdi['country_fe']})
res9 = reg9.fit(cov_type='cluster', cov_kwds={'groups': high_hdi['country_fe']})

# Regression models, very_high_hdi
reg10 = smf.ols(formula='change_co2 ~ change_ln_gdp + C(country_FE)', data=very_high_hdi)
reg11 = smf.ols(formula='change_co2 ~ change_ln_gdp + ghg + C(country_FE)', data=very_high_hdi)
reg12 = smf.ols(formula='change_co2 ~ change_ln_gdp + ghg + change_ei + C(country_FE)', data=very_high_hdi)

# Fit the regression models
res10 = reg10.fit(cov_type='cluster', cov_kwds={'groups': very_high_hdi['country_fe']})
res11 = reg11.fit(cov_type='cluster', cov_kwds={'groups': very_high_hdi['country_fe']})
res12 = reg12.fit(cov_type='cluster', cov_kwds={'groups': very_high_hdi['country_fe']})

# Create Stargazer object
stargazer = Stargazer([res, res2, res3, res4, res5, res6, res7, res8, res9, res10, res11, res12])

# print(stargazer.render_latex())

# Load image of regression results table without Country Fixed Effects coefficients 
Image(filename='reg_table2.png')

Image(filename='reg_table3.png')

from sklearn import metrics
from sklearn import tree
from sklearn.ensemble import BaggingClassifier, RandomForestClassifier, BaggingRegressor, RandomForestRegressor, GradientBoostingRegressor
from sklearn.metrics import mean_squared_error,confusion_matrix, classification_report

x = reg_df[['change_ln_gdp', 'change_ei', 'change_hdi', 'ghg', 'hdi_cat', 'country_fe']]
y = reg_df['change_co2']

sqft_tree = tree.DecisionTreeRegressor(max_depth=3).fit(x,y)

y_pred_tree = sqft_tree.predict(x)
print('Mean Squared Error Machine Learning:', metrics.mean_squared_error(y, y_pred_tree))
# Print Mean square error of res3
print('Mean Squared Error Preferred regression:', ares3.mse_resid)

Mean Squared Error Machine Learning: 0.2661688150766863
Mean Squared Error Preferred regression: 0.2716194813783243

sqrf_fig = plt.figure(figsize=(25,20))
sqrf_fig = tree.plot_tree(sqft_tree, feature_names=x.columns, filled=True)

# define and fit
regr_RF = RandomForestRegressor(max_features=5, random_state=1).fit(x, y)

#predict
pred = regr_RF.predict(x)

#calculate MSE
mean_squared_error(y, pred)
Importance = pd.DataFrame({'Importance':regr_RF.feature_importances_*100}, index=x.columns)
Importance.sort_values('Importance', axis=0, ascending=True).plot(kind='barh', color='r', )
plt.xlabel('Variable Importance')
plt.gca().legend_ = None

	Country	ISO	Year	Total CO2	CO2 Per Capita	GDP Per Capita	HDI	Energy Intensity
0	Afghanistan	AFG	1990	2.024326	0.189281	NaN	0.273	2.429634
1	Afghanistan	AFG	1991	1.914301	0.178155	NaN	0.279	1.154009
2	Afghanistan	AFG	1992	1.482054	0.122916	NaN	0.287	0.723151
3	Afghanistan	AFG	1993	1.486943	0.106182	NaN	0.297	0.925850
4	Afghanistan	AFG	1994	1.453829	0.094065	NaN	0.292	1.113139

	Country	ISO	Year	Total CO2	CO2 Per Capita	GDP Per Capita	HDI	HDI Category	Energy Intensity	ID
0	Albania	ALB	1990	5.520602	1.675415	2549.746801	0.647	1	2.599493	ALB1990
1	Albania	ALB	1991	4.290226	1.299249	1909.319160	0.629	1	3.010896	ALB1991
2	Albania	ALB	1992	2.517116	0.761899	1823.503609	0.614	1	2.821280	ALB1992
3	Albania	ALB	1993	2.337999	0.708332	2057.692048	0.617	1	2.245501	ALB1993
4	Albania	ALB	1994	1.923806	0.584034	2290.143917	0.624	1	1.941011	ALB1994

	Total CO2	CO2 Per Capita	GDP Per Capita	HDI	HDI Category	Energy Intensity	GHG Policy Ratified
count	3520.000000	3520.000000	3520.000000	3520.000000	3520.000000	3190.000000	3520.000000
mean	245.356864	4.652458	15305.409897	0.672432	1.488636	1.533110	0.305114
std	891.076734	5.490882	17067.118842	0.174236	1.163471	1.418874	0.460521
min	0.000000	0.000000	287.661076	0.216000	0.000000	0.209657	0.000000
25%	4.138031	0.618075	3131.659673	0.532000	0.000000	0.833408	0.000000
50%	25.465533	2.368603	8808.760671	0.694000	1.000000	1.168157	0.000000
75%	147.809827	7.290669	22244.178287	0.813000	3.000000	1.836268	1.000000
max	11472.369171	33.296139	131630.774358	0.962000	3.000000	23.142410	1.000000

	rank	country	cumulative emissions per population 2021	rank_0	country_0	cumulative emissions per capita 2021
0	1	Canada	1,751	1	New Zealand	5,764
1	2	United States	1,547	2	Canada	4,772
2	3	Estonia	1,394	3	Australia	4,013
3	4	Australia	1,388	4	United States	3,820
4	5	Trinidad and Tobago	1,187	5	Argentina	3,382

Greenbacks and Carbon Tracks: the Evolution of Economic Development and CO2 Emissions¶

Project 1¶

Introduction¶

Variables¶

Data Cleaning/Loading¶

Summary Statistics Tables¶

Plots, Histograms, and Figures¶

Project 2¶

THE MESSAGE¶

Maps and Interpreations¶

Project 3¶

Potential Data to Scrape¶

Potential Challenges¶

Scraping Data from a Website¶

Scraped Data Visualizations¶

Adding a New Dataset¶

Final Project¶

OLS Regression Models¶

Machine Learning¶

Conclusion¶

Citations¶

Datasets¶

Literature Review¶